function y = eq_func(x,Nb,pi_m,func_var,opt)
%===============================
%
%===============================
beta=func_var.beta;i_m = func_var.i_m;
alpha1=func_var.alpha1;alpha2=func_var.alpha2;alpha3=func_var.alpha3;
lambda_func=func_var.lambda_func;dlambda_func=func_var.dlambda_func;
rq_ratio=func_var.rq_ratio;cost=func_var.cost;df=func_var.df;

if opt==1
    y(1)=alpha1*lambda_func(x(1))+alpha3*lambda_func(x(2)+x(1))-i_m;
    
    psi = beta*(alpha2*lambda_func(x(2))+alpha3*lambda_func(x(2)+x(1)))+beta;
    dzdD=-alpha3*dlambda_func(x(2)+x(1))/(alpha3*dlambda_func(x(2)+x(1))+alpha1*dlambda_func(x(1)));
    dpsi = beta*(alpha2*dlambda_func(x(2))+alpha3*dlambda_func(x(2)+x(1))*(1+dzdD));
    rho=x(3);
    xi = (1-rq_ratio)*max([1+rho,1/(1+pi_m)])+rq_ratio/(1+pi_m)-cost;
    y(2)=xi*(dpsi*x(2)/Nb+psi)-1;
    y(3) = df(psi*x(2)*(1-rq_ratio),Nb)-1-rho;
elseif opt==0
    y(1)=alpha1*lambda_func(x(1))+alpha3*lambda_func(x(2)+x(1))-i_m;
    
    psi = beta*(alpha2*lambda_func(x(2))+alpha3*lambda_func(x(2)+x(1)))+beta;
    dzdD=-alpha3*dlambda_func(x(2)+x(1))/(alpha3*dlambda_func(x(2)+x(1))+alpha1*dlambda_func(x(1)));
    dpsi = beta*(alpha2*dlambda_func(x(2))+alpha3*dlambda_func(x(2)+x(1))*(1+dzdD));
    rho=x(3);
    xi = (1-rq_ratio)*max([1+rho,1/(1+pi_m)])+rq_ratio/(1+pi_m)-cost;
    y(2)=xi*(dpsi*x(2)/Nb+psi)-1;
elseif opt==3
    psi = beta*(alpha2*lambda_func(x(2)))+beta;
    
    dpsi = beta*(alpha2*dlambda_func(x(2)));
    rho=x(3);
    xi = (1-rq_ratio)*max([1+rho,1/(1+pi_m)])+rq_ratio/(1+pi_m)-cost;
    y(1)=xi*(dpsi*x(2)/Nb+psi)-1;
    
    
end